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^ : Abstract 

^ - The saturation of the tearing mode instabihty is described within the standard framework of reduced magne- 

■^ohydrodynamics (RMHD) in the case of an r-dependent or of a uniform resistivity profile. Using the technique 

o f matched asymptotic expansions, where the perturbation parameter is the island width the problem can be 

r^iplved in two ways: with the so-called flux coordinate method, which is based on the fact that the current profile 
I ' 

^s a flux function, and with a new perturbative method that does not use this property. The latter is applicable to 
C^ore general situations where an external forcing or a sheared velocity profile are involved. The calculation pro- 

Q-i' 

^ides a new relationship between the saturated island width and the A' stability parameter that involves a.\nw/wo 
• ^erm, where wq is a nonlinear scaling length that was missing in previous work. It also yields the modification of 



^l}e equilibrium magnetic flux function. 
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In magnetized plasmas with current inhomogeneities, magnetic islands can develop as the consequence 
of tearing mode instabilities. This mode corresponds to a global magnetic perturbation that is resonant 
on a magnetic surface where its wave-number is perpendicular to the magnetic field. From a theoretical 
viewpoint, the instability is allowed by non-ideal plasma effects such as resistivitj*^, electron inertiai^ and 
more generally kinetic effects^"^. On the resonance surface, there is a modulated current sheet whose 
current density is proportional to the classical tearing mode stability parameter A'. This parameter, 
whose precise definition is reviewed below, depends solely on the ideal magnetohydrodynamics (MHD) 
properties of the system. Instability occurs when A' exceeds a critical value, (A' > Ac > 0). This value 
depends on the actual physics included in the model. For instance, in visco-resistive MHD, Ac is a function 
of viscosity and resistivity^. Above the critical value, the resonant magnetic surface can break up (tear) 
and is substituted by a chain of topologically distinct structures called magnetic islands. Until recently, 
tokamak operation avoided the formation of such islands, since they result in greater radial transport 
and hence deteriorate particle and energy confinement. However, in certain experimental conditions, a 
magnetic island may help the plasma form a stable internal transport barriep^. In the reversed-field pinch 
(RFP), the occurrence of several magnetic island chains leads to magnetic chaos with reduced confinement, 
but the formation of a single chain is desirable, since it should provide good magnetic flux surfaces and 
a laminar dynamo. Therefore, a correct description and understanding of the nonlinear tearing mode is 
both important for thermonuclear fusion and for advancing the theory of plasma self-organization. 

In the theory of magnetic island formation, nonlinear effects come into play as soon as the island width 
exceeds the width of the boundary layer at the resonant surface, as given by linear theory. The nonlinear 
tearing mode is classically described by applying resistive reduced magnetohydrodynamics^ (RRMHD) to 
the model of a static plasma slab, in the limit of small dissipation. The magnetic island region is considered 
as a boundary layer whose nonlinear features are dealt with, while the outer region is adequately described 
by linear theory only. The inner and outer solutions are then matched asymptotically. Rutherford^ showed 
that the island growth is sufficiently slow that inertia can be neglected in the inner solution. The fiuid 
equations of RRMHD then reduce to a mere force balance (Grad-Shafranov) equation. Rutherford also 
showed that the island width w grows with a linear time dependence in the early (small island) nonlinear 
phase. A quasilinear calculation predicted a further nonlinear slowing down of this growth^. 

The saturation of the tearing mode is a difficult issue, and its solution has been a stepwise process 
covering almost three decades. In 1977, a seminal work using a quasilinear calculation provided a first 
version of the formula A'{wsat) linking A' to the saturation amplitude Wsat of the island widtb^i. In 1981, 
a new technique was introduced to deal with the case of a non vanishing current gradient J'^qirs) on the 
resonant surface, and gave the correct expression for the leading order term in Wsat of /S.'{wsat) in this 
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term and finite [3 corrections were also given in Ref. llJ. In 2004, the rigorous expression for /S.'{wsat) 
was provided for the case J'(,g{rs) = (Harris sheet pinch) by two different techniques, one of which uses 
explicitly the fact that the current profile is a flux function^, while the other does uoIf^. 

The aim of this paper is to provide a deeper insight in the nonlinear tearing mode by introducing, in 
particular, a powerful perturbation technique which does not use the fact that the current proflle is a flux 
function. This makes possible, whenever useful, to avoid the assumption that plasma inertia, pressure or 
viscous effects are negligible, which opens a new route to deal with a background velocity proflle and/or 
an external rotating forcing. Furthermore, whenever the current proflle is a flux function, we show that 



the technique introduced in Ref. can be simplifled to make the derivation more direct, which divides 
by more than two the necessary algebra. As in all the above quoted works, Wsat is assumed to be a small 
parameter in both of our techniques, and A'{wsat) is computed in a perturbative way for the case of a 
cylindrical geometry for any current gradient J'^^{0). 

Reference 113 provides a short presentation of the new perturbative technique in slab geometry by the 
same authors. Our flnal formula for A'{wsat) depends on Wgat through a term Wsat^^{'UJsat/u)o) where wq 
is a nonlinear scale length which was absent in previous work. We show that the nonlinear tearing mode 
comes with a modiflcation of the background magnetic flux. It should also come with a modiflcation of the 
resistivity proflle. Indeed, since temperature is typically uniform on a given magnetic flux surface, so is 
the resistivity. Therefore, the development of a magnetic island is bound to flatten the resistivity proflle. 
Since our RRMHD model does not incorporate the evolution of temperature, we provide the formulas for 
the saturation of the nonlinear tearing mode for both uniform and non uniform resistivity proflles. The 
same formula for A'{wsat) was obtained in parallel by another group who applied the technique of Ref. 
12 to higher orde r^^'^'^ ; refill also provides expressions for A'{wsat) where the feedback of the magnetic 
island on the resistivity proflle is accounted for in the case of a time independent thermal conductivity. 

This paper is organized as follows. In Sec. m we introduce the main equations and we specify the 
different normalizations and notations adopted throughout the paper. The solutions of the (linearized) 
outer equations are given in Sec. IIIII We then tackle the problem of the saturated tearing mode in the 
nonlinear inner boundary layer (Sec. llVjl . where the flux coordinate method, which exploits the fact that 
the current is a flux function is described flrst (Sec. IIV Bjl . We outline how the calculations can be carried 
out consistently to the leading signiflcant order in the island width expansion. This provides a benchmark 
for the more general perturbative approach, the main focus of this paper, which is presented in Sec. Fvl 
Section ED yields the modifled Rutherford equation providing the whole nonlinear evolution of the island 
width, as well as the (lowest order) modiflcation of the equilibrium magnetic flux function; it also provides 
a brief discussion of the validity limits of our method. Finally, Sec. IVIII is devoted to the conclusions. 
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A. RMHD equations 



The RMHD equations are given by: 



dtA±ip + (v.V) A^^ = B.VJ + uA^Lf 



(1) 



dtiJ + B.V<^ = r]{Jeg - J) 



(2) 



J 



(3) 



where the mass density is uniform and taken equal to 1 for simplicity, ± denotes the plane perpendicular 
to ez, ip is the poloidal flux (i.e. B = B^ez + V x (^ez)), (p is the stream function (i.e. v = x Vip), v is 
the viscosity and r] is the resistivity. In the following, two models are considered: model A, in which the 
resistivity is uniform [r] = tia), and model B, in which it is not, but the electric field is {rjB{r)Jeq{r) = Ez). 
The equilibrium current profile Jeqi^^) then fully determines the plasma equilibrium which we have further 
assumed to be static {ipeq = 0). 

B. Single helicity perturbation of the equilibrium 

In the remainder of this paper, we work in cylindrical geometry and for now only consider the saturation 
of the single helicity (m, n) perturbation of the equilibrium, to reintroduce time dependence at the very 
end of our calculations, which, as we shall see, can be done very simply. Therefore, all quantities henceforth 
depend on two variables only: r, and r = mO — nz/R, where R is the (simulated) major radius. Equations 
(IH) and ^ then take the following form: 



where q = rB^/RBg is the safety factor, and, for any functions / and g, [f,g] = drfdrQ — drgdrf is their 
2D Jacobian or Poisson bracket. 




(4) 




(5) 
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It is easily seen that 'B.Vi/j ^ 0, and, therefore, ip is not appropriate to describe magnetic surfaces. In 
contrast, the hehcal flux function ip* = ip + BzTir'^ /2mR verifles B.V^^* = 0, since B = Bzh + V x (^*ez), 
where h = + nr/mR eg. We thus from now on work with ip* instead of ip. Finally, introducing the 
following normalizations 

r = rof ; e = ^ ; J = Jo J ; ip* = /iorQ JoV^* ; ^ = ] V = VoV, (6) 

where ro is the minor radius, ^o^oJo = B^, rjo = rjA for model A and rjo = E^/Jo for model B, equations 
(jH) and © along with Ampere's law can be written in the following way: 



777 

-[^,r]=V{Je,-J) (8) 



J=- (-drirdrr) + ^d^rr] + 2" (9) 

where we have omitted the "~" for the sake of clarity. S = v^rofiQ/riQ is the Lundquist number. Re = vaTq/v 
the Reynolds number, and va = Jo'To^fflo the Alfven speed. Since S ^ 1 and Re ^ 1, equation (0) merely 
gives^: 

[J,^*] = (10) 

This means that J = Ji^ip*), which is the basis of the flux coordinate method. Equations ((HI), © and (|T(Hl 
are the basis of the following analytical work. 

D. Boundary conditions 

Equations © and (fTOll must be complemented with boundary conditions. Cylindrical geometry 
requires: 

lim (^^^ = drr\r=o = = 0) = 0) and Jim (^^^ = dM^^^ = (v^(r = 0) = 0) (11) 

Furthermore, since the plasma boundary is taken to be at r = 1, the normal velocity component should 
vanish on this surface: 

dMr=i = (12) 
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to carry out the following calculations. Suitable boundary conditions are, for example, drip*\r=i = 0, as 
in Ref. or drip*\r=oo = 0, as in Ref. 0. 

This ends the definition of the differential problem to be solved. 

III. OUTER SOLUTION 
A. Outer equation 

As usual in tearing mode theory, we solve the differential problem as a boundary layer problem by 
matching an outer and an inner solution, and we approximate the outer solution by the linear ideal one. 
We set: ip* = ip*{r) + S'^ipipir) cos(r) + o(5^). Linearizing (fTTHl yields: 



From now on, we assume that m > 2, since we use the "constant-?/^" approximation which does not work 
for the m = 1 mode. Equation (fT3|) then has two regular singular points: one at r = and the other 
at r = r^, where ipe'qi'i"s) = 0, or equivalently qeqi^s) = fn/n, since Qeq = rt/{nre/m — ip*'g). is the 
location of the {m,n) rational surface and, in the frame of the "constant-?/'" approximation, we can choose 
i^iFi^s) = 1 SO that 6 is the square root of the perturbation amplitude at r = r^. The indicial equation for 
both points shows that, whenever the boundary conditions (fTI|l are satisfied, ipip behaves like r™ around 
r = 0, and has a logarithmic singularity at r = r^. This means that our perturbation expansion breaks 
down near the rational surface, which eventually has to be resolved thanks to a boundary layer centered 
upon it that we henceforth refer to as the "inner" region. 

B. Inner limit of the outer solution 

Let p = r — Vs, we now give an expression for ipip close to p = 0. To do so, we first expand Jeq, ip^q 
and qeq as Jeq = Ylixd'^iP^ i^eq = Ylii>2^ip\ ^^'^ leq = Ylii>o^iP\ whcrc the a's, h' s and c' s are of course 
related to each other. In particular: 




(13) 



ne ao 



m 




(14) 
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S' ± A' r S' ± A' / 1 



\ ( rr? 



a 



+a<;p+-(l 



pnin|p|+0(p3) 



(15) 



where ± = sign{p\ a = — ai/262, and /? = — (ai/462?"s + 02/^2 + 0^/862)- A' and S' are two constants 
that are determined by the boundary conditions. They can be expressed as: 

A-=hm f^'^Firs + e)-.^' ir.-e)\ = lim f ^^^^^^ + " - 2.(1 + lne)l 

(16) 

A' is the usual tearing mode stability parameteA Note that, contrary to A', S' depends on normalization. 
Indeed, had we chosen to normalize r with respect to, say, instead of tq, S' would have been changed 
into S' = S' + 2aln|r^/ro|. This remark will prove to be important when we come to the saturation 
equation. 

The logarithmic term appearing in (fT^ implies that the perturbation expansion breaks down in a region 
centered on the rational surface. This comes from the fact that the quasilinear term in (|TT11l becomes of 
the same order as the linear one when p is sufficiently small. Indeed, it is easy to see that: 



[lpeq,S JifCOSt] + [6 IpiF COS T, Jeq] = 0{6 ) and [6 IpiFCOSTjS Jifcost] = 0{ 2, 



(17) 



where Jii^cosr = — A(^ii;^ cos r), and Eq. (fTTjl immediately shows that the quasilinear term is no longer 
negligible when p ~ 5. Therefore, in order to deal with the boundary layer, we use the stretched variable 
^ = p/6. Since the outer solution 'iplgir) + 6'^ipiF{r) cos(r) + o{6^) is going to be matched to the inner 
solution which will be computed in the ^ variable, we now re- write it in terms of ^: 

'S'± A' 



h^" + a 



S'± A' 
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+ ^(l- — )(ln|e|+ln5)e^cosr 
2 V ar 



&3r + 
1 

ar 

■2 



2 
1 

2'^ 



+ a(ln|^| +\n6) ] ^cosr 
1 

2 ~ ar, 



1-— ) + ^(— + /3-^(3-— )) le'cosr 



(18) 



where q = sign(h2). Note the appearance of In 5 terms that are due to the logarithmic singularity. 



IV. INNER SOLUTION 
A. Inner equations 

Until otherwise stated, we now consider model A (constant resistivity) and will address model B later, 
where only minor modifications will occur. We work in r) variables. Matching with the outer solution 
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5^ < 1, Eqs. © and © become: 

f 1 - 5^] [C, + o{5') = J-{ao + ai5^ + as^'^') + o{5^) (19) 



= ^ - ^slc + r f 1 - 4) + S *'9.'C [ + o(i=) (20) 



where the Poisson bracket is now taken with respect to (^, r) variables. 

As regards the inertia equation, some care is needed before simply re-writing equation (fTfl|l instead 
of 0. Indeed, we first have to compare the order of magnitude, in the boundary layer, of the different 
terms appearing in the latter equation before making any simplification. To do so, we need to have some 
information on ip. This can very simply be done by writing ip = d'^ipipir) sinr + o(5^) and solving ^ to 
order 5^: 

Tj [m %l)[p riplp 



ip*^ \r m m J 

We see that, since (pip goes as p~^, then, in the boundary layer (p ~ 5), = 0(1). It is now straightforward 
to show that the terms on the right hand side of (0) can be neglected if the following conditions hold: 

5 > S~^'^ and 5 > {S.Re)-^'^ (22) 

The first of these is basically that 5 be greater than the resistive layer width, a fact already pointed out in 
Ref. 



and the second one is that it be larger than the visco-resistive length. Provided (f22)) is satisfied, 
the inertia equation can be written, in the inner domain, as in equation (fTIHl : 

[C,J]=0 (23) 

In the following, we introduce two independent calculations of the nonlinear inner solution, which 
eventually has to be matched to (fH^ . using a perturbation expansion in 5. The first one basically replaces 
equation (f^ with J{C,,t) = j(C(^)T)), i-e. it uses the fact that the current profile is a fiux function. We 
only give a brief account of it, since it is an improvement of the technique described in Ref. Q, and gives 
the same solutions as the second one which is new. The latter is indeed based on a classical perturbation 
expansion of all the involved functions, and is more fiexible since it works also for problems where J is 
not a function of ip. 

B. Flux coordinate method 

As already mentioned, equation implies that J(^, r) = j{C{C,, r)). This makes natural the following 
change of variables: 

(e,r)-(C,r) (24) 
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result, a given value of C may correspond to several values of ^ for a given value of r (at the maximum order 
of our calculations, there are only two values corresponding to the topology of a classical magnetic island). 
Therefore, the change of variable TIM is one to one only in local domains. In such domains, we can solve 
Ci.ii r) = constant for ^, which yields ^ = r). Consequently, ip{^, r) becomes r) = r), r). 

It follows that (fT9|l can be recast into: 

^1l[d^x(l + 5-\\ dMC.r) = m - Je,{5X{C,r)) (25) 



The periodicity in r of $(C, t) implies that the integration of (f25|l on a flux surface covering the — vr < r < tt 
interval, which we refer to as S'^, gives : 



A closed flux surface is described by a series of functions X{(,t) (two for the maximum order of our 
calculations) . These functions enable the generalization of Eq. (f^ to closed flux surfaces by interpreting 
the integrals as loop integrals on a given flux surface. 

Equation (f26|l is the fundamental equation of the flux coordinate method. It was already derived in Ref. 



l^for the case of a given resistivity proflle (i.e. model B). The inner solution can be derived by combining 
it with a single differential equation, Eq. (pUl) . through the following iterative procedure. The first step 
sets 5 = in Eq. (|^ . This yields jo = «o, which is set in Eq. ^Jl^ to provide Coi = |^2|^^ + + &(t) 

where a(r) and 6(r) are two unknown functions which are determined by matching with the outer solution 
(fT8|l . Since the next order in the expansion of C provided by Eqs. (|26|l and (f2n|l is 6, the matching brings 
to a(r) and 6(r) terms of order 1 and 6\n6. As a result, Co* = Co + 'iaSlnSC, cosr, where 

Co = \b2\e + ^cosT. (27) 

Then the leading orders Xoi of X can be computed by solving ( = Coj(^oj(C? t), t) at orders 1 and 5 In 5 



for Xqi. This yields Xoi(C,T) = Xq — qa6ln6cosT, where Xq = ±,/\h\^ — (jcosr). This calculation 
requires |Xoi(C,T)| ^ 51nl/(5, which excludes a small neighborhood of ^ = 0. 

This ends the first iteration of the calculation. We notice that two values of ^ are related to one value 
of C at this level of approximation. The next iteration starts by setting Xoi in Eq. (f26| . which brings 
orders S and 5^ In 5 to j. These orders are brought into Eq. (pill) . which brings contributions of orders 
6 and 6^ In 5 to C which are completely defined by matching with the outer solution, and so on. The 
exclusion of a small neighborhood of ,^ = is required at all orders. The calculation brings an expansion 
of A' in 6 which is provided in the next section. We notice that each step of the fiux coordinate method 
approximates 5^ by its expression given by the available approximation of X{(,t). 
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results of the flux coordinate method provide perturbation expansions ( = J2i^''Cii i(C) = 



Jo = flo 



ji = ai Xodt-Xo dr 



32 



(9^Xo dr 



^ { («2 + X^d<:Xo - ji (^d^X, + d^Xoy-^ + ai9c(XoXi)| dr 



and 



df^Xo dr . 



(28) 



Xi = -9cXoCi(Xo(C,r),r) 



(29) 



The expression for is given in the next section. 

Notice that a differential equation for X can be obtained by rewriting Ampere's law (f2n|l with respect 
to the new variables: 



2ne 



6 1 - 6X/r, j 1 



drX 



dr 



drX 

Wx 



(30) 



The flux coordinate method brings a series of simplifications to that of Ref. (i) no Ansatz is made 
about the solution, which brings only the non vanishing orders in 6; (ii) the use of the ( variable enables 
the same calculation to be formally done for flux surfaces inside and outside the magnetic island, and 
simplifles the calculation of flux surface averages; this divides the necessary algebra by more than a factor 
two; (iii) fewer quantities need to be deflned to proceed with the calculation. 



V. PERTURBATIVE METHOD 



A. Zeroth order 



The method is simply based on equations (fT9|) . (f20|) and (f23|) along with the following natural ex- 
pansions: C = Co + ^Ci + ^^C2 + o(5^), J = Jo + SJi + 6^J2 + o(5^) and ip = (po + 6(pi + o{6), where a 
quantitiy with index n has an order smaller than 1/6 and larger or equal to 1. From (fTTI|l . we directly 
obtain Jq = ao, which, using (f^ and matching with (fTH| . immediately yields (f?7| . We recognize the 
"constant-?/'" approximation that is valid to lowest order. 

In the rest of the calculation, it will prove most useful to work in (Co, t ; ±) variables where ± tells 
the sign of ^. Let /(^,t) be any function of the old variables, then we should introduce /(Co,t;±) 
such that /(Co(^, t), T ; ±) = /(^,t). Nonetheless, to simplify formulas, we do not make that distinction 
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e = ±v/|&2|-nCo-^cosr). 

Finally, we define contours as: 

ifx>l C± = {(e,r)GR±x [-7r,7r]/Co(e,r)=x} ^^^^ 
if - 1 < s < 1 = {(^,r) G M X [^^^ + arccosx, -arccosx]/Co(^,r) = a;} 

where the first line describes open curves and the second closed ones (from now on, we omit the ± 
superscript for open curves, which should not make any confusion). These contours merely represent 
lowest order magnetic surfaces. In the following, we systematically make use of the fact that, for any 
single- valued and r-periodic function /(Coi ^ ; ±): 

drfdr={drf) = Q (32) 
B. First order 

From (f23|l . it can readily be seen that Ji = ji(Co ! =t)- We note, already, that we use the same notation 
ik as the one introduced in Sec. IIVB| since it will be shown shortly that they do refer to the same 
functions. Writing (fT9|l to order 5: 

Tfi 

-(,—2\h2\idripQ = ai^ - ji (33) 
rs 

and integrating the equation above along Cq,, we immediately derive the following expression for ji: 

J,{Co;±)=2na^HiCo-l)/{r') (34) 

where H is the Heaviside function. We now make an important remark concerning ji. Indeed, expression 
(|3i|l has a derivative singularity at Co = 1? which is not physically acceptable. As already mentioned in 



Ref. ll^, this problem can be resolved thanks to a thin boundary layer centered around the separatrix 
Co(^;X) = 1- However, contrary to what is claimed in Ref. [3, it is not inertia but viscosity that is 
no longer negligible in equation (f^ (see Appendix A). The current profile could thus, in principle, be 
regularized by solving the problem in this secondary visco-resistive boundary layer, using once again the 
technique of matched asymptotic expansions, a procedure similar to that already performed in Ref. [21I to 
regularize Rutherford's solution^. Nevertheless, that treatment needs not be done explicitly for present 
purposes. We therefore assume that ji is regular from now on, although we always use expression ()34|1 in 
our calculations, which makes sense as long as it does not lead to divergences. Finally, using together 
with (pn| provides (i. 

Ci = 7^( r n{x ;±)dx-^ r ^r^^ dx] - + Air)^ + i?(r) (35) 



2^2 \Ji ' Ji ^(x,r;±) / 3rs 
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Before proceeding to the matching with (fT8|l . we note that, given any function f{C,o{^,T),T]±), it is 
possible to obtain its asymptotic expansions in one of two ways: either express it as an explicit function 
of t) and directly derive its expansion as |^| ^ 1 (i.e. expand the related /(^,t) function), or expand 
it as a function of (Co, ^ ; ±) for Co ^ 1 and, only then, re-write that expansion with respect to (^, r) 
while making |^| ^ 1 (see Appendix B). We use the latter method since it is much more convenient to 
implement. 

It is straightforward to prove the following expansions: 

±1 / / COST 

ji(x;±) = ^^v^+0(x-^/^) 



^1 ic^ x3/2 



= (2Co) 



Jl dx = ^^ (2Cof^ ±Q + 0(C(r^/^) where n= Urn ( ± f^" ji dx - , 

Jl 3v/2|6^ Co-- V Jl 

j ^ = aiCo + ^cti In -\/Co cos r — ai + S(r) + 0(C(7^) where S = j — ai — qai—^ — ^ dx 

(36) 

where VL = ujai/ a/2|62|, and is a numerical coefficient that is approximately equal to —1.54. 
Then, taking (pIHl) and (jT^Kl) . we simply set Co = 1^21^^ + ^cosr and expand for ^ 1: 

-I (M) - - - >" i« i) « — « - ± ^ - -^(^'^ - - <'^) 

Since, in our approach, we have only taken into account the zeroth and first Fourier components from the 
outset, we should match only these two in our calculations. Thus, making use of ((H]), we can match (irTTIl 
with ((IHl), which determines A{t) and B{r) 

A(r) = a-^ f "^'^^^^^ ~^ -aln(5 ) COST and 5(r) = (38) 

\ 2 2o2 / 

where Si = tt^^ J^^S(r) cosrdr. What is important to note is that the first A' term appearing in 
(fTR|l cannot be matched with Ci- Indeed, it would require the inclusion, in A(r), of a quantity of the 
form ±^A'cosr/2, which is not allowed since it would lead to Ci's being singular at the rational surface. 
Therefore, A' has to be matched with higher order terms, which will precisely provide the saturation 
condition we are looking for. Note also that, for the same reason, we cannot compensate the ±$7/262 term 
with B{t). Since it is not matched in ()18p either, it implies a modification of the equilibrium magnetic 
fiux that is of order 6^, which will be discussed later. 
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Moving on to the next order, it is easy to show that (f23|) now implies J2 = j[{(o ', ±)Ci + j2(Co ; =!=)• Here 
again, j2 has to be determined through equation (fTOll . To order 6'^, we have: 

-2 



TTl I ^ \ 



Besides, (f33|) gives: 



-<i^2\h2\dr'^Q = cti - 7- and - q—2\b2\dc_adr'^o 



Ji 



and, using (|in|) . (f39|) can be recast into: 



2|&2|e e 



Ci 



J2 



-ai<9coCi+ ( — + ^2 1 ^ 



Then, integrating (jlT|l along gives the expression for j2 



32 



Ci 



(39) 



(40) 



(41) 



(42) 



Note that j2 turns out not to depend on ±, i.e. is even in ^. We also see that, to know J2, we needed 
the full expression for (^1. However, since we stop the calculation at order 5^, only (9^^2 is now required to 
later perform the matching. We therefore integrate (piUl only once: 



216^ 



'21.3 



3r2 



JiCi C: 



— + ^cosr + — / 



Ci 



- ^ ) rfx + C(r) 



(43) 



2^2'^ Ts n ^^2 ^fcosrx 

where, again, C(r) has yet to be determined. 

Now that we know both the inner and outer solutions, we can proceed to the matching procedure, 
which is actually the most difficult part of the calculation. We begin by expanding the first terms in (|i3|l . 
which is easy to do since we have already derived (f37|l : 



(32 



^<^^^ + y + 6^ + i262 



1 



y + a(ln + In 5) V K cos r 



(262)^ 



^cosr ±a— 1 +C r + — 

rl 262 V CirJ 262 



Co 



? cos T 



3idx 



Ci 



— ^ 1 



(44) 



The last term of that expression is, therefore, the main part of the calculation. 

We first determine its diverging part. Making use of (f34| . (j35|l and (|42| . it is possible to derive the 
following expansions: 



■ r. (C,i\ Tai ( ai l\ r— a\ 



at cos r 



(2I&2I) 



-3/2 



aif2 
4&2C0 



+ 0(Co"'/^) 



J2 
1 



262 



"2 + — + ^ ) 1^ - ^^(2|62|)-^/^ + 0(Co-^) 
^ 1 dx 



g c os T 

1&2 



Ci 



(45) 



at 



_^ 02 

1262 3 6rJ I62I 



Co 



7a? a2 



ai 



2462 ' 262 ' 1262r. 



cos T > + convergmg term 
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ai 



462 rs 



- — -a 1 — + a(ln|e| +ln5) ^cosr 



±a-?- fl - — V ± lim 



Co 



0,1 do 0-1 

— ^ + — H 

1262 3 6r, 



Co 



+ 



Co 

f cos T 

2 



Ci 



Tar 02 



Lieu lei 



— -rrr 1 fix 



2462 262 1262^^ 



cos r 



0(1) 



(46) 



Using (fni and the definition of /3, it is possible to show that the first line of this equation matches all 
diverging terms of (fT8|l . 

The last expression on the right hand side of (|i6|l is the "converging term" mentioned in (|^ . In order 
for the matching to be complete, we have to select its first two Fourier components. We begin by the 
zeroth Fourier component 



± lim 



Co^oo I 47162 



We first note that 



Co 

dT \ dx \ jidx 

' <; cos T 



± 



Co 

dr I dx [ jidx 

' f cos T 



Ci 



32 



Co 



a; 



]_ I ^ _|_ '^1 

126^ 3" 6^ 



Ji 



dr I dx { h aidxCi 



Lieu lei 

where use has been made of and We then derive the expansions below: 



ai 

02 H 



3/2" 



(47) 



(48) 



-1 
47r62 
1 

47r62r 
ai 

47r69 



cir / dx^ = =F';^ f-f|^ 

<;cosr -J VP2I 

Co 



3/2 



+ 0(11 



dr 



dri i,,=±,|L(|L]''^±JL + „,l 

fcosr -J^s VP2I/ ^02 '^s 

^« / ai a2 \ / A„ \ 3/2 

(ix 9^Ci = T'i ( ^ 1- 

\Qrs 1262 



(49) 



I&2I 



+ a- + o(l) 



and thus see that the zeroth Fourier component exactly cancels out the ±a(l — l/ars)fi/262 term on the 
second line of The final matching condition is therefore given by C(r) = and 



A' = 6 lim , 

Co ^00 I 7r02 



Co 

drcosT I dx ( jidx 

' ? cos T 



lel 



J2_ 

lel 



2 ( + — + ^1 
2462 262 1262r, 



A. 

I^2| 



(50) 



The limit on the right hand side can be evaluated numerically, which eventually gives: 

a2 , fliS' ^ a2(ln ln5) \ 1.39 0.65 ai 



A' c^dl 1.64 



+ 



o{S) (51) 



&2V2I62I 462v/2|6^ (2|62|)5/2 J bl,/2\h\ r,h^/2\h\ 

This equation is the final result of the matching procedure and will be re-written in a more pleasant 
form shortly. We see that (f50| requires the knowledge of divergences so as to substract them from the 
current profile integral. Furthermore, checking that the diverging terms do match correctly with the outer 
solution (fT8|l gives confidence in the validity of the saturation condition (f5T|l . Lastly, given (f28|) . (f29|) and 
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expansion of (f28|) yields the same current profile functions Jn that were derived above. In particular, the 
jfc's introduced in Sec. IIV Bl and the ones derived in Sec. El are indeed the same functions, which justifies 
the use of the same notations, as was said previously. Therefore, the flux coordinate and perturbative 
methods do provide the same result, as already claimed. 



D. Non uniform resistivity 

We now treat model B (non uniform resistivity) using the perturbative method (it is easy to derive 
the equivalent of (f26| for that model) , and show that only slight differences occur while solving the inner 
equations. Of these, only Ohm's law is changed into 

<,!^5 (l _ si) [C, 4>] + o{6') = vmJ- 1 (52) 

where ri{p) = J2i>od'iP^ satisfies rjJeq = 1. Therefore, the a's and d's are easily related to each other and, 
in particular, do = a^^, di = —ai/a^ and d2 = aQ^(af/ao — 02)- 

The only difference then appears at order 5^, where j2 now becomes: 

-^^((^--^«^-(^^" 

Using (jHHll . it is straightforward to prove that has the same asymptotic expansion as that given in (jlHl) 
for j2- Consequently, the only differences that we have to determine are those coming from equations ()47p 
and (f5n|) . Making use of (|^ . one easily shows that: 

lei; 3 6rJ\\h\J J 262 V^s " aoj 

(54) 

Therefore, 8^(2 now has an asymptotic term of the form ^aiQ/2b2ao which we cannot compensate with 
C(r). It actually has to be matched with the order 6^ correction to the equilibrium magnetic flux that we 
mentioned earlier and which we will deal with somewhat later. It is interesting to point out that such a 
term was forbidden in model A because of current conservation during the relaxation (see Appendix C), 
and our calculations are thus consistent. 

As to the flrst Fourier harmonic, it can be shown to merely include a new term in ()51|) whose coefficient 
can again be computed numerically and is approximately equal to —0.28 (5a^/ao&2 a/2|62|- 



Ci 



- — + — )- ai^^oCi 



'in 



(53) 



lim 

Co ^00 



47r62 



dr 



(,0 
dx 

C COS T 
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A. Saturation equation 

The matching conditions that we have obtained for both models actually provide saturation equations 
for the island width, which is defined as w = 46/ a/2|62| (i-e. it is the width of the separatrix of zeroth 
order magnetic surfaces as described by Co)- We first need to reintroduce time dependence so as to make 
comparisons with previous results easier. Equation (fT9|) is modified by adding —dt{6(t)^Q to the left hand 
side, whose lowest order term is 255 cos r, which changes (j3i|l into 

Ji(Co,t;±)= (27raii7(Co-l)-2ao5(|^^)/(r'> (55) 

Since the new part added to ji is even in ^, it contributes to the matching with A'. We therefore take 
this new term into account and rewrite (f5T|l with respect to w and for both models: 

aow = 1.22 A' + w\-a (a\nw + —] - 2.21 + 0.40 — - ^ - 0.17 Aa—1 + o(u;) (56) 
1 2 V 2 y r, 262 ao J 

where A = for model A and A = 1 for model B. We now reintroduce normalizations explicitly and show 
that 



^ ^ J'egjrs) / _ 2 



B _ a2 _ J':,{rs) A 2 



2 2rlh2 2J,,(r,) V' s 
d fiQ 2 d 

dt Veqirs) at 



(57) 



where s = r sq'f,q{r s) / qeq{r s) is the shear parameter, and use was made of the fact that 262 = 0-05/(2 — s), 
which can be proved using relations (fTlj) . Equation along with (fHTjl . then provides the final evolution 
equation 

^"-1.22A' + ^.(:^('^ln^ + ?V2.21^^ + 0.40:^ + f +0.17A^Uo(^) (58) 



VeqiTs) dt ' I2V ro 2 J ' ' 2 ' 2-s 

where the A' and S' are now to be understood as dimensional parameters whose definitions can trivially 
be deduced from (fT6| . A result similar to (f58|) was obtained using a variant of Thyagaraja's technique^^, 
the only difference being the numerical coefficient of the wA/vs term which was there found to be approx- 
imately equal to 0.22. 

We now make an important remark concerning ([Hl^ . Indeed, it can easily be shown that, despite the 
Inw/ro term, it does not depend on the normalization length tq. The reason for that is straightforward, 
and comes from the comment made below (fT6|) : S' precisely depends on tq in such a way that the 
combination A\n.w /tq + S'/2 is actually normalization independent. It is therefore natural to define an 
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yj^^roexp-^- where S' = lim f ^^^^^^^±^^±4^^^^ - 2^(1 + In -)1 (59) 
and eventually rewrite (fHI^ in an explicit normalization independent way: 

It is interesting to review past work in the light of this consideration. First, one observes that the results 

n n n 

obtained in Refs. HU and HJ do depend on length normalization. The work done in Ref. is formally 



normalization independent, since a change of normalization length induces a change in the saturation 
equation of order S, higher than the order 5 In 5 at which the calculations were stopped. However, the 
normalization issue strikes back in any practical application of the formula when an explicit value of 6 must 
be set in the logarithm. Finally, the equation given in Ref. Q turns out to be normalization independent, 
but lacks important order 6 terms and does not provide any result for the symmetric case (i.e. ^ = 0) 
either. 

We finally mention that all the calculations which have been carried out so far can very easily be 
modified so as to fit the slab geometry case already treated in Ref. [3 One simply has to apply the 
following simplifications: 

771 

ao = 1 , 62 = — 1/2, ^ 00 , s— >-oo, ifj* ^ if) , p^x and r^x (61) 



In particular, (jHUll immediately allows us to recover the evolution equation derived in Ref. 







B. Modification of the equilibrium magnetic flux 

As pointed out when we derived the asymptotic behavior of the inner solution, there has to be an order 
modification of the equilibrium magnetic fiux in the outer solution for the matching procedure to be 

complete, which is a result different from that obtained in Ref. Q, where it was found to be of order 6^. 

We therefore write ip* = iplq{r) + 5'^7pip{r) cos(r) + 5'^7pQp{r-) + o(5^). From what we have done earlier, we 

know the following conditions that TpQp must satisfy: 

lim ^pir) = ±7TTrT and lim tp'opir) = ipA-^— (62) 

r^Ts ^1^2 1 r—>rs ^|02| CLq 

We now have to derive an equation for ipop and we will see that our simple perturbation technique is 
much easier to implement than the one used in Ref. Q which is based on the J = J{ip) property. Setting 
(p = 5^v9ii?(r) sinr + 6^ipQF{r) + o{6^), we see that equation (|TT1|l is trivially satisfied to order 6^ and thus 
write ((HI) at that same order, which merely gives Tp'^p + Tp'^p/r = and, consequently: 

Tpopir) = D + E\nr (63) 
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violate the condition described in Appendix C. Physically, it does not lead to any change of the equilibrium 
magnetic field but merely reflects the fact that the saturation of the magnetic island leads to an increase 
of the (normalized) poloidal flux per unit length that, to lowest order, is equal to 26^Q/2\b2\- Making use 
of the normalizations (jH]) and integrating on the whole cylinder finally gives the total change of poloidal 
flux below: 

-OM8i2oRJ'e,{rs)wlt (64) 

where Wsat is the saturated island width. 

In the case of model B, things are not so straightforward. Indeed, in order for (j63|l to match conditions 
!2jl . one would naively write 

,5 , Q airs, r 



The problem is that such a solution would not satisfy the boundary conditions of the problem given in (|TI| 
and would actually be singular at r = 0. This means that the second condition in (|62|) is not acceptable. 
That problem can be solved by setting C(r) = aiQ/2b2ao in which is perfectly allowed. ipQp then 
flnally becomes 



which results in the poloidal flux being changed into 

-0.048 /.oi?j;(r.)^L f 1 - '-^j^ In ^"j . (67) 

V •Jeqys) '^s J 

This change of order in the flux is a natural consequence of Jq = ao and of the expression (f3i|) for 
ji. Indeed, they tell that, in the nonlinear regime, the current has a plateau inside the island. Since the 
width of the island is 0(5), this brings a change 0(Jgq(rs)5^) to the magnetic field and 0( Jgq(rs)5'^) to 
the magnetic flux. 



C. Validity limits of the method 

Since we have used a perturbation expansion in 5, the first condition that should be met for our 
calculations to be valid is obviously 5^1. For instance, in the case when ^ = and -B > 0, equation 
(|60| would predict exponential growth of the island width. However, since our result is no longer valid 
when b approaches unity, this does not necessarily mean that the system would lead to a disruption, for 
there might be saturation with a large island. 
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away. Suppose that A' is such that A'5/2 ~ 1. Then, in (fT8|) . the lowest order term should include the 
A'S one and, therefore, Co would become 

Co = |&2|e' + ^fl±^eVosr (68) 



2 

However, this would not be allowed, because of the singularity at r = r^, and our method would thus lead 
to a dead end. 

This basic analysis is in good agreement with recent numerical results obtained by Loureiro et alr^ for 
the symmetric tearing mode in slab geometry. Indeed, since, in that case, the island width is simply given 
by w = 45, the condition A'(5/2 ~ 1 gives w/S! ~ 8, which is surprisingly close to the condition derived in 



Ref. [2^ (wA' ~ 8.2). When this condition is met, Loureiro et al. observe the formation of current sheets, 
which, typically, means that the first Fourier harmonic of the perturbation is no longer dominant with 
respect to higher ones, which contradicts the "constant-?/^" approximation used in our approach. 



VII. CONCLUSION 



We have provided a rigorous solution to the simple tearing mode problem in cylindrical geometry using 
both the fiux coordinate method and a new perturbation technique, and our calculations can directly be 
transposed to the case of a plasma slab. The final evolution equation contains all terms of order w and 
has been explicitly shown to be normalization independent, a necessary physical requisite. We have also 
shown that the saturation of the tearing mode leads to a modification of the equilibrium magnetic fiux 
function which we have been able to fully determine and which is consistent with the condition of current 
fiux conservation in the case of a uniform resistivity profile. Lastly, we have discussed the limits of validity 
of our approach and we have derived a qualitative condition on A' which is in good agreement with the 
recent numerical study carried out in Ref. [2^ 

Besides the actual results shown here, it is important to appreciate the importance of establishing solid 
analytic techniques for future work. In particular, the perturbation method, which does not rely on a 
functional dependence between the current and the fiux function, is rather promising for the treatment of 
more general models than conventional reduced MHD. For instance, the forced tearing mode in rotating 
plasmas is currently being revisited and a first application to the static case was given in Ref. Q. In the 
longer term, two-fiuid models with diamagnetic effects should also be accessible to analytic investigations 
with these techniques. 
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APPENDIX A: VISCOUS BOUNDARY LAYER AROUND THE SEPARATRIX 

We want to evaluate the first neglected terms in (f23| and derive their behavior around the separatrix 
Co = 1. Equation (f33|l gives: 

Furthermore, using ([2T| . the matching condition for is 

avg sinr ,^ 9, . ^ . 

which is automatically satisfied by the first term on the right hand side of (jXlJ. 
To lowest order, the corrections to (f23|) are given by: 

-^5^[Co, Ji] ~ ^[<^o, dl^o] - ^J']^^^A ^o (A3) 

We thus see that, if the magnetic Prandtl number S/Re is of order unity, the main correction is due to 
viscosity and we can neglect inertia. Integrating ()A3|) on C^^ then gives the following condition: 

si (eolfo) = di m dlt,) + {(i' [ dl iy)) = (A4) 

W . . e.. . _ (e[.i {1} - I ™ .... ...e < - 0. ...0. 

neglected term in (OTHl is given by: 

4/)2 /r \ 2 / r ( 11 



S.ReS^ \mJ \ Jo ^° U 
It can be shown that, close to the separatrix, ji behaves as 
27ra 



Ji 



U\n 32 - In (Co - 1))"^ + 0(Co - 1) (A6) 



2 



If the position is not too close to the 0-point or to the X-point, it is then easy to see that 



20 



5„ 



3/2 



Co-l-[jj (A8) 

where 5^ = {S.Re)~^^^ is the visco-resistive length. (IA8|1 then determines the size of the viscous boundary 
layer centered on the separatrix which should, in principle, regularize all the singularities that appear in 
the current profile. 

APPENDIX B: VALIDITY OF THE Co ^ oo ASYMPTOTIC MATCHING 

We first define the equivalence relation below: 

Let G (M«)', f^^.g^ lim (/(x) - (^(x)) = (Bl) 

a;— >±oo 

If we define the function /(Co,t;±) such that /(Co(^, t), t ; ±) = /(^,t), what we want to show is the 
following property: 

If / and / then ^o(Co(^, r),T]±)^ f^{^, r) (B2) 



Let e G M"*"*, then, by definition of the equivalence relation (|BHl . we have: 

3M,GM+, V(Co,r)GM^ Co>M,^|/(Co,r;±)-/oo(Co,r;±)|<£ (B3) 



Now, let = \J\h2\-^{Me + 1), then it is clear that: 

V(e, r) G R\ \i\ >M,^ |/(Co(e, r), r ; ±) - /oo(Co(e, r), r ; ±)| < £ (B4) 
Since /(Co('C) t), t ; ±) = f{^, t), we immediately deduce (IB2jl . 



APPENDIX C: CURRENT FLUX CONSERVATION IN THE CASE OF UNIFORM RESIS- 
TIVITY 

In model A, Ohm's law is written as: 

777 

-b,V'*]= Je,-J (CI) 

Integrating the left hand side on a poloidal section 5 = {(r, r) G [0, 1] x [0, 2n]} gives: 

— [(p,ilj*]dx dy = — dr dr {dr{ip* dr4>) — dr{ip* dr(l))} 
s ^ Jo Jo 



dr {r{0,T)d^(P\r=o - r{l,r)dr(^\r=i} + [ dr{r{r,2n)d,4\r=2.-r{r,0)dr(^\^=o} (C2) 
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boundary conditions (fT2|) . the last line of ()C2j) is equal to zero and, therefore, integrating the right hand 
side of ()CHl gives: 

(Je, -J)dxdy= [ (Beq - B) .rfl = - f / dra,V'1r=i - 27r<;(r = 1)) = (C3) 
where we have used Stokes' theorem and again the periodicity in r. Condition (l('3|l thus imposes that 

dr^OF\r=l = 0. 
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